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ABST RACT 


The mathematics of extrapolating known statistics of 
components to the probability density function of a system's 
performance measure 1S considered. Quadrature sum inte- 
gration schemes for evaluating the resulting required inte- 
gration are examined, and alternate integral approximation 
schemes are developed utilizing Monte Carlo methods. A 
Simple electrical circuit example illustrates the use of 


these techniques. 
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in ERO DUCTION 


The desirability of being able to statistically predict 
a system's performance measure 1S obvious when one considers 
the mass production of systems. Since performance measures 
of certain values may not be desirable, it is generally un- 
profitable to produce such systems when the probability that 
these undesirable values can exist is beyond the limits 
determined by the requirements of the situation. 

Presently performance data for systems are obtained in 
several ways: (1) A worst case analysis [l] to determine if 
undesirable performance will occur within the parameter toler- 
ance; (2) A moment method [1] in which the performance 
measure is sampled so that its mean, variance, and higher 
order moments are determined for approximation by known dis- 
MiaPowemon funetiGas, (3) A functional formulas method [2] in 
which performance measure's distributions are obtained from 
breaking performance measures into a series of products and 
sums for individual integration; (4) Monte Carlo methods 
Pyelo teri whlch performance measure distributions are gener- 
ated by creation of random processes by which the parameters 
of the random process lead to an approximation to the desired 
pels tee ut On 

The purpose of this thesis is to investigate the mathe- 
Matics required in obtaining a performance measure and at- 


tempt to arrive at a reasonably fast and accurate method of 





extrapolating the statistics of a system's component variables 
to create the distribution of the desired performance measure. 

ene II examines the mathematical form and statistical 
tee OF the performance measure for implicit and explicit 
functions and presents two digital computer peacraenchern 
schemes uSing Gauss quadrature sums for evaluating the re- 
Seeing required integration. 

Chapter III investigates crude Monte Carlo methods for 
determining performance measures, and the confidence limits 
involved in the method. Then Monte Carlo schemes to estimate 
Single integrals are examined. The multiple integral exten- 
fens of the single integral estimates are then derived. 

Chapter IV discusses the relative merits of the methods 
outlined in Chapters II and III, and gives a comparison of 
the results of these methods applied to a simple example con- 
Eioaeing Of an electrical circuit in which the component 


volMesware the random variables. 












Pie otal tel iCAL ANALYSIS 


The probability density function of a system's perform- 
ance Measure can be determined mathematically from knowledge 
of the functional relation between the peeeserereinice measure 
and the system component-variables, and knowledge of the 
Merle probability density function of the component- 
variables. Most often this requires integrating a complex 
function over several variables. Although the integration 
can be performed analytically in some cases, digital computer 
integration schemes must normally be utilized. 

This chapter examines the mathematics of determining the 
Dmobablility density function of a performance measure, first 
for performance measures which are explicit functions of the 
component variables, then for implicit performance measures. 
Two integration schemes are then presented for performing the 


mieeg@ation On A digital computer. 


Pee exe LECLIT CASE 

Consider a performance meaSure, Z as a function of the 
system's component-variables, ee where x" is the transpose 
of an n-vector of elements Xa Xos Sees xX: It is assumed 
that the X'S are statistically independent, so that the 
omen eebablivty density function of es is given by the 
Peeqdmermoer tne probability density functions of the individ- 
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Wal xX. 's. 
L 





ieee beobaemtity density function of Z = g(x") can be 


written as [7] 


E Th 
p,(z) = \, yarns 72) AY (1) 


fk 
where Bye oS By Pama 7OlmceorOobabilaty density function 


for X and Z, V is the n-dimensional space determined by the 


component-variables, g is the functional relation between Z 


and the component variables, xe anGdwedVr——dx.dxk. ...dx 


eZ 


mirc ansrortmacwon Of varltables 1s defined such that 


a 4 | en 
Y (X5,Xa,¢---7%)) > xX a lox 


ah 
x aa) (2) 


T 


Xy OZ *) (3) 


where f£f is the functional relation between the component- 
Varlable xX) and the performance measure Z, equation (1) 


becomes [7, 9] 


py(z) = [ pytlt(a.y™) y"1[o(2) [au (4) 
U aes 


where U is the n-l dimensional space determined by the n-l 
: d's _ a 

Cemeenent Variables Y , dU = dy,dy5---dy,_j5 Ys = Keays 

Meroe lacobian Of Ehe transformation (2) and (3). 
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Tome canmeomwritten as 








: OX, OX) OX) 
ce ey ¥n-1 
IX. 
OZ . 
dl (44) 
Ox ox 
_ 
OZ oy ea 
Ox 
al 
oa 0 0 
a Ox 
ms oe 1 
= ; Ce (5) 
Ox, 
OZ 


Wo@eresteis the n-l by n-l identity matrix. 

Siiees ese JOInt probability density function of statis- 
tically independent random variables is the product of the 
meeryvrqual probability density functions, equation (4) 


becomes 
py(z) = | py Celery") sy") [oxy/d2] py? (y"au (6) 
a el = 


In the case where the performance measure 1S a known or 
explicit function of the component-variables, the inverse 
Functional relation f relating xy to Z can usually be ex- 
pressed analytically as can the partial derivative of Z with 
respect to xX) Gu cdiveor cwem@tier Component-=varlables. Cir- 


cumstances may require reordering the component-variables 


MOG om@cmunmmenesce relatiensceim a convenient form. 


Lat 





ecm este yecOlsStacr the Series R-L—-C circuit in 


megure sw. The current I through the circuit is 


= Ei) 
(jw) = FFL —17ae7 a 


and 1S a maximum for wL = 1l/wC. The angular frequency at 


Pmeen the current falls to 3 db of the maximum occurs when 
Ree lie —- 1/WC (8) 


Oi, 


_ Boy 
pena eee 2 BR /4i> + 1/LC (9) 


Let the performance measure Z be the upper 3 db frequency 
and xt = (R,L,C), so that 
PP x) XK. -/4X.° + 1/xX2Xx (10) 
ot al 2 a 2 23 


amaertrom (8), 


_ Ty = 
oe 04 a 2X5 = 12x, Gi) 


so that the partial derivative is 


/92 SG yee: 


5 (a) 


5 


If each of the components has a Gaussian distribution 
about its mean such that 
exp [-(x.-x.)°/20.°] 
By 2 Sayer ae (13) 
5 (27 ) Oe. 
J 
where x, is the expected or average value of the variable 


Bac and o. is the standard deviation of the variable X.. 


aE 





L C 


T0000 i. 
Lbs w) 
(SY Ej) R 


Figure 1. Series R-L-C Circuit 
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Substituting for Xy Pilate e— |, Py [f(z,y )+y) 
ib 
becomes 


exp [- (zx, -1/zx4-%,) °/20 4°] 
Py 2 i (G4) 
dl (ae, oT 


Since X, and X., are independent, pyt(y) tween mtOduCL OL 


w 5 
Py (x5) and Py (x4), and (6) becomes 
2 3 
— ,2 — ,2 a2 
~(zx,-1/2x,-%)) (x5 -X,) (x, x) 
exp [—— >. - — > 
2 2 2 
oF po 20 20 26 
a 
ZL i 3/2 5 a 
ee, (27) 04050311/(x5+1/2-x,)] 
(15) 


With a-priori knowledge of the probability density func- 
tions of the component-variables, (6) or (15) can be inte- 
grated on a digital computer using one of many algorithms 
avallable. The two quadrature formulas presented below offer 
routines which are fast, and for the number of points at 
which the integrand is evaluated, two of the most accurate 


mewitrnes avVallable [5, 8]. 


Pea LePEtecry CASE 

If the performance measure iS not known explicitly, but 
Pm—enewarmmea Mmpltcrtly from a solution of a mathematical 
model of the system, relation (3) cannot be determined. One 
may, however, resort to the cumulative Clcotri DUE mon Oocarmned 


Peon ly) by [7] 
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ee 





Z 
ea es < 2) = | | PyT,, (x ,z) avVdz (16) 
2 A te 


— OO 


Since 
De (x 12) = Bee 272 pT (x") (17) 


where Po xt (2/x) Poeeneweconaltzlonal probaesz lity of Z given 


that a Naswocecumued; and for a given a Py xT (2/x) is 
elther 0 or 1, the cumulative distribution can be approxi- 


mated by 
AL 
F(z) = Y Az | DyT(x ) av (18) 
Zee V-=- 
where the integration over the space V may be performed 


utilizing the quadrature formulas below with the value of 
T Ak E ; 
DyT (x Peoomcalcnlatcad for Z2(x )< 2 and PyT (x ye == 7 Org ak te 


Z(x') > z. 

The cumulative distribution function for the performance 
measure can then be fitted to a polynomial, and if desired, 
the polynomial can be differentiated analytically to obtain 


EieompboOoabllity density function. 


C. QUADRATURE SUMS FOR EVALUATING INTEGRALS 

The form of the performance measure has been reduced to 
eMeemeeagral or mMultaple integral which, in general, cannot be 
evaluated in closed form. All algorithms which approximate 
the value of an integral by a linear combination of values of 


the integrand are exact [5] for the integrand being of a 


15 





certain degree or less; and in most cases the degree is one 
less then the number of points at which the integrand is 
evaluated. Quadrature sum algorithms utilize orthoganal 
polynomials to arrive at the form of the linear combination 
of values of the integrand; and as a result Breweaee for the 
integrand being a polynomial of degree one less than twice 
the number of points at which the integrand is evaluated. 

[5] 


Consider the integral 
b 
= | fetes) gi) ax (19) 


where a and b are real numbers, finite or infinite, g(x) 1s 
PEE atetrary funNct1On and h(x) 1S a particular weighting 
function to be described. Then by selecting a polynomial 
Q, (x) of degree n such that A(x) Q, Cx) 1s orthoganal to xn 
moe ali m= 0, 1, -.., n-l, so that 


b 
| xh (x)Q_ (x) dx = (0 (20) 
a 


Lf (20) is approximated by a linear combination of the 
integrand 
b rm, n - 
| xe (x) Q) Cx) dx = d Ay Xy QO, Cy) (21) 
a k=1 
then the right hand side of (21) will be equal to zero for 
any set of values of A. 1f the X'S are chosen such that they 


are the zeros of Q, ts). 
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By forming the n sums 


b n 
| x h(x) dx = ) ee 22) 
a k=L1 

Mommie = O, b,...,n-l, and for the n values of Xp being the n 


zeros of Q(x), the n values of AL. can be found such that 
b n 
| h(x)g(x)dx = ) Ag (x) (23) 
a k=1 
is exact for all polynomials of degree 2n-1. [5] 
It can be shown that if the polynomial QO, (x) 1s known 


then the values of AL Saneoe found to be [5, 8]. 
= 2/(1-x 2) [ao (x \ fax] 7 O24) 
“a k n ‘*, 


1. Gauss Legendre Quadrature 


Micm=mcsemObVYLOUS Chol1ce Of w(x) 18S 1, and if the 


2 


feeerval fa, b] can be normalized to the interval [-l, 1] by 
a change of variables, (23) becomes 
1 1a 
| g(x)dx = d Ag (xy) C25) 
“51 k=1 


and the polynomial Q, (s) is a Legendre polynomial of degree n 


[Se 


QA 2 n 
Q, (x) = d=) (26) 


dx™2™"n1 
Ihe n values of xX, can be found from the zeros of 
Q,(s) and the values A. from (24) or the values can be ob- 
PieneGderconelaole Ili for values of n from 2 to 6 and from 
lbule tor evaloues Of nm £rom 2 to 48. The magnitude of the error 


bound can be shown to be [5], 
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2n 
e es (me max| 24 


 Onalnempe x dx" (27) 


2. Gauss Hermite Quadrature 
As 1S Often the case when dealing with probability 
density functions, a variable iS GausSian in distribution, 
and as a result an integral of the form 


roa) at 
| h (x) exp (-x*) ax = | AA Cxy, ) (28) 
k=1 


must be evaluated. This is the form of (19) where w(x) = 
exp (-x*) . 
The class of Chebychev-Hermite polynomials are 


orthogonal with respect to exp (-x*) Over the real line, and 


are defined by [5] 
n 2 q” 2 
Q, (x) = (-l1)” exp(x”)—T[exp(-x”) ] (2) 
n 
6b 
The values of x, can bewweOunameLoOmerEne Zero’ Ss Of 
(29) and those of AL from (24). However these values have 
been tabulated and are included in Table I for values of n 
mae nte Go and in [5] for values of n from 2 to 32. 
The magnitude of the error bound for approximating 
(28) by a Gauss Hermite Quadrature sum can be shown to be 
fl 
1/2 max) d°"h (x) 
on 


Ss dx 


ae, 


= | Gere 
Zee! 


BeeoeDplitecation to the Example 
The integrand for the example, equation (15) is of 


the Gauss Hermite form Since the variables all have Gaussian 
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TABLE I 


GAUSS HERMITE QUADRATURE 
COPY ECIENTS AND 4EROS 


Number 
of Evaluation LEOLOSe x: Coefficients A. 
Points, n J J 


2 +0.707 1068 0.886 2269 
3 +1.224 745 0.295 4090 
0.0 1.181 635 

4 +1.650 680 0.081 31284 
40.524 6476 0.804 9141 

5 +2.020 183 0.019 95324 
¥0.958 5725 0.393 6193 

6 +2.350 605 0.004 530010 
+1.335 849 0.157 0673 
+0.436 0774 0.724 6296 
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TABLE IT 


GAUSS LEGENDRE QUADRATURE 
COEFFICIENTS AND ZEROS 


20 


Number 
Of Evaluation Zeros xX. Coefficients A. 

Pores, n J 
Z Oe! 7 omy 3 1.000 0000 
3 uO. 7/4 596 7 O25 55175 556 
0.0 0.888 8889 
4 ple SOL eG 3 0.347 8548 
BO. 339. 9810 Oe65Z2. 1452 
5 a0. 90a 7 98 0.236 9269 
Oso 204693 0.478 6287 
0.0 0.568 8889 
6 =029 32 4695 Only) 3245 
FO-o6) 2094 0.360 7616 
t Opec O92 0.467 9139 





GtStyPoueions. By changing variables, let 














' Qe 2 
Xo = (31) 
420 
2 
X.—-X ; 
sy 
x, = (32) 
/20 
S 
: dx. 
aX, = (33) 
Y20 
Z 
dx 
dx. = (34) 
V¥20. 


Since X, and X. are independent, the double integral 


becomes the product of the two integrals and (15) takes the 


ete 
| | g(2,x,,x3) exp ( Xo -X3 )Axj dx 
n n 
=e AA G(Zp%> 1X3) (35) 
j=1 k=1 J 25 3k 
where 
g(2,X5,X3) 
. 1/2, = 1 te Z 
7 expi-t2(2 o., X5+K, yz 2 O4x4+x,) X13 /20 5 ] 
1/2 3772 il 25 ae ria -1 
2 (1) Oy [2 J5 x5 +X, 2) fae (2 O3X4txX3) | 
(36) 


The value of n can then be chosen as desired and the sum 


in (35) formed using the values of AL and x and x from 
a1 3x 


Table II or [5]. The expression (36) becomes quite formidable 


Pak 


[ee 





to differentiate after two or more differentiations, so that 
one must assume that the error bound (30) is small. The 
oxi mation (35) has been performed for n = 6 and the mean 
Values R = 1000 ohms, L = 1.0 millihenry, C = 0.001 micro- 
farad. The standard deviation of each variable was assumed 
to be 3.33% of its mean value. The results are contained in 
G@rapter IV. 

Not all variables have Gaussian distributions, and most 
performance measures, 1f functions of many variables, are 
implicit functions. If the example performance measure is 
considered to have been obtained implicitly, with the same 
density functions for the components, the performance meas- 
ure's cumulative distribution function can be determined by 
approximating (19) by a Gauss Legendre Quadrature sum. Since 
wot =+.5) 1S very nearly zero, the limits of -m and © can be 


changed to -30 to 30 for each variable and by the following 





change of variables, normalized to the interval [-l, l}. let 
x, -K; 
t — — 
x, = ane vee lp oS (37) 
1 
t — 


so that (18) becomes 
ii 1 
Eee) = ) | | Azg(x),X5,%3) dx) dx5dx, (39) 
Z 
neh a ee | 


where 


Ze 





"Zz a ve _ a 
27exp(-x) iD Xo UD Xs Te) 


g(x) ,X4,x4) =  — > > (40) 
aay (27) 3/2 


The approximation of (39) by the summation 


n Lp gl 
B{zZ) = ) ) ) ) ee (KS kh ea XL) Ca} 
2<2 i=l j=l k=l = + pre FUT, 182573, 


just requires the n values of A. and the n values of xs and 
k 


evaluating g at these values of Xe -. It must be remembered 
k 


that g takes the form of (40) only for K(Ce | z for some 2; 
@enerwise g(x') = 0. 

Equation (41) is evaluated for n = 6 and the same values 
of the components as above for the Gauss Hermite example, and 
the results are. contained in Chapter IV for comparison. As in 


the Gauss Hermite example, the error is assumed to be small 


hero thne Lactorials in the denominator of (27). 
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PEE. MONTE CARLO METHODS 


This chapter examines the most crude of the Monte Carlo 
techniques for obtaining statistics on a performance measure, 
to methods which attempt to approximate single integrals in 
Such a manner that the approximation iS reasonably accurate. 
The approximations are then extended to apply to multiple 
integrals in such a manner that fewer evaluations of the 
integrand are made than in conventional methods. 

Monte Carlo methods consist of solving problems of a 
computational nature by constructing a random process for the 
problem, such that the parameters of the random process are 
the desired quantities of the problem. 

Random processes usually imply random variables which 
must be generated in a manner as to represent typical values 
from the variable's probability density function. Random 
number generation is then an important aspect in any Monte 
Carlo method of approximation. Although random number gener- 
ation is a field of its own, Appendix A treats the subject 


Suitably for purposes of this thesis. 


Poo ce RUDE MONTE CARLO 

iiemuesteaeraide  tornm of Monte Carlo is that of taking sam- 
ples of the performance measure to obtain an expected value 
and perhaps a variance of the performance measure. If the 
component-variables are generated in such a way so that the 


values are representative of their distribution functions, and 
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for each set of values we of component variables generated, 
the performance measure Zs = Z(X7) 1s evaluated, then the 


expected or average value of the performance measure is [7] 


Ho = E(Z;) = (42) 


oF Te 
te et F) 
i) 


where n is the number of samples. 

Ditto Tomesutticiently large m, Z iS an accurate estimate 
of the average value of the performance measure, no knowledge 
Smeene form Of the probability density function of 2 is de- 
rived; and consequently no knowledge of the likelihood of 
other values of Z 1S obtained. 

A more useful scheme would be to order the sample values 
of Zs as follows to obtain a cumulative distribution for the 
performance measure. 

Consider a random variable S such that if a value of the 


performance measure Z 1s sampled, and the sample value Z. 1s 


less than an arbitrary value, Say a, then Ss = 1+ and 1f£ Z. 
1s greater than the value a, Ss = 0. Then S represents a 
Eandom process. If the values of the component~variables ea 


are sampled from their density function as above, then the 


eumubative distribution function for Z can be approximated by 


F(a) = P[Z<aJ = 


B= 


n 
va (43) 
where n is the number of times the performance measure is 
Sampled. 

Since each of the component-variables are sampled inde- 
pendently for each sample value Ze the number of samples is 


not directly dependent on the number of variables. 
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This crude method yields a cumulative distribution func- 
Mmwonewilch may be fitted to a polynomial and differentiated 
feo yrically to Obtain the probability density function for 
the performance measure Z. 

As an example of the two methods, consider the example 
Gaeenapter Il. Since the three variables X 


X51 and X 3 are 


a 


assumed to be Gaussian, three variables Wir Wor and W, are 


generated by the method in Appendix A to approximate Gaussian 


Sister ibutions, each with O mean and variance of 1. The val- 
ues of the samples Xa Kos x. are then 
x, = WO, + ree a Ole ear: (44) 


where O; and Xs are the standard deviation and expected val- 
ue, respectively of the variable Wi. These values are then 
used to evaluate the performance measure. If the expected or 
average value of the performance meaSure is desired, the cal- 
culated values of the performance are averaged as in (42). 

If a cumulative distribution is desired, then for k values of 
a j=1,2,...,k the random variable S defined above is summed 
Zo. £Orm PEL) c These values for the examole are listed in 


@aapteer IVY, ana a computer program listing for the computa- 


EWOmeis provided in Appendix B. 


lou CONFIDENCE INTERVALS 
The crude Monte Carlo method above, as with all Monte 
Carlo methods, yields results which are as accurate as de- 


Sired within a probability or confidence determined by the 


desired accuracy, and the method of approximation used. What 
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ZHOLlOws 1S a Statistical estimate of the error introduced by 
use of Monte Carlo approximations. 
If Ais an event (such as the event that the values of a 


performance measure is less than some specified value), and R 


is a random variable, such that R, —i Pee eOocecurs On the ith 
sample, and R, = 0 otherwise, and if T is a random variable 
such that 
n 
a (45) 
i=l 


for n samples, then T is the number of times that the event A 


Occurs inn samples. The expected value of T is 
sk) = me (R.) = np (46) 


where p is the probability of the event A occurring. The 


waeiance of T is 
Tek) — 1. (AG) 


where 5° is the variance of the event A. 


From Chebychev's inequality [7], 
2Z 2 
P = P[|T-np|< ad] > 1 - o”/nd (48) 


where dis an arbitrarily small number. Equation (48) states 
that for a fixed confidence of |T-np| being less than a small 
value d, the number of samples to achieve the accuracy, d, 
1s bounded by 

Z 


ee (49) 


BS lD) 


where P is the desired confidence or probability. 
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Another way of looking at (48) is that for a fixed con- 
fidence P, the error of the estimation varies as on Ye tt 
the variance of the random variable can be reduced to 0, the 
confidence of the estimation could be 1 for an error one (tl, 

For crude Monte Carlo, the random variable S, which takes 
Smeeclues Of | whenever the sampled performance measure is | 
less then a specified value, and values of 0 otherwise, has a 
binomial distribution. If p is the expected value of S, then 
the variance of S is se = p(l-p) [7]. Thus for a confidence 


eee 0.9 Of being within d=0.01 of np for crude Monte Carlo, 


Oe te wesc a p(1-p) 10° (50) 


Pemece thie maximum value of p{l-p) is 0.25, inequality (50) 
states that n must be greater than twenty-five thousand in 
@ieer tO insure the accuracy of 0.01 with a confidence of 


ie. 


C. VARIANCE REDUCTION BY ESTIMATING INTEGRALS 

Romvocmocenmin Chapter £, the difficulty of determining a 
performance measure's probability density function lies in 
evaluating the integral (1). Two routines were presented 
which approximated the integral with some accuracy. The 
crude Monte Carlo method illustrated above can be thought of 
as an approximation to an integral. Other methods exist, 
however which reduce the variance by hundreds of thousands 
over crude Monte Carlo [3]. These methods are presented below 
for single integrals and then a general extension is made for 


multiple integrals. 


23 


_— 





1 
fits the. estimator for the integral, | f(x) dx, where 
0 


XxX 1S a variable normalized to the interval [0, 1], then by 
letting 


ee) 
if = =— ee 


where Py (y) 1s the probability density function of the ran- 


dom variable Y and f(y) is the value of the function f at 


some point y sampled from Py ly), the expected value of I is 


p(t) = eff) = [ syyyay (52) 
Py (y) 10 : 
[fiers then necessary to select the probability density func- 
EON Py (y) in such a manner as to minimize the variance. 
Define a function f£*(x) such that the function which is 
fexy-t*(x) 1S a straight line passing through the end- 
Pemmes Of £{(X) On the interval [0, 1], as shown in Figure 2. 
That is, 
£¥*(x) = f(x) - (1-x)f(0) - xf(1) (S38) 
The difference function f(x) -f£*(x) can be readily inte- 
grated and can be used as a first approximation to 
1 
| ies rc as 
0 


ut 
| (ene ee ecm nl 72 £(0) 4 1/2£ (1) (54) 
0 


Now if I* is defined as 


[* = £* (y) (55) 
Py Cy) 





it becomes necessary to sample f*(y) to get the estimate from 


(soi ts4)eand (55) £or LI 
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£(X) 


(x) 
% 
(X) 
Figure 2. Difference Curve Estimation. 
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eee eee (iy) £(0) = y(2) 
a9 — rr eee 
5 : Dy (y) 7 (56) 


To determine Py ly), Terex) be quadratic in x, that is 


Zid) & pneo ee eure | (57) 


Peemoubestituting for £ in (56) the £f of (57), and replacing I 


fee Ehe integral of f in (57), that is 
—T = A/3 + B/2+C (58) 
one can obtain for Py ly). 


Py(y) = 6y(l-y) (59) 


It can be shown by substitution that the approximation I in 
(oa), With Py (y) Seno et cmonzerOo Variance Estimator for 
the integral of quadratic integrands. 


The function £*(x) could have been defined as 
ae lV) a ave (0) =7(l=y) £( 1) (60) 


with a resulting estimator of 


+ S a) eas ee . Seley (O}e— (l=y)£(1) 
2 2 6y(l-y) 


(61) 
WimemmcmaylSO exact fOr £(x) quadratic im xX. Since y and 
l-y have the same distribution, but lie on the opposite side 
of 1/2, the average of (56) and (61), as can be shown by 


pubstimutrion, Fesults in a zero variance estimator for £(x) 


Sore in x{/3|]. Thus 
eee Oy tf (1) f{y) + £(l-y) - £(0) - £(1) 
i ) 7 'S RS as Oe) 


I, meme waot for £ (xX) seu bilc ville. 


od. 





In the same manner as above, by defining f* as 


ie Ore ee £ (x) ~ (1-x) (1-2x) £(0) +x (1-2x) £ (4) -4x(1-x) £(5) 
(63) 


a difference curve which is quadratic and passes through the 
two end points and the center of f(x) becomes the basis for 
the initial estimate for the integral. Applying the methods 
above, a zero variance estimator for quartic integrands be- 


comes [3]. 


a £(0)+£(1)+4£(1/2) 
4 Ba..6, 
he img Claes eo) ty i-2y yk 4y (1-y) £(1/2) 
30y (1-y) (1-2y)7 
(6 4) 
and a zero variance quintic estimator becomes 2 
: £(0)+£(1)+4£(1/2) 
5 6 
, £(y)+£ (-y) - (1-2y) * [£(0) +£ 1) J ~ey (1-y) £172) 
30y (1-y) (1-2y) * 
(65) 


iiemerebaotiity density function for both the ela Gdieaiz ie 
pmmGleeclbic zero Variance estimators is (59) while that for the 


Slsere and quintic zero Variance estimators is 


Py (y) = 30y (1-y) (1-2y) 7 (66) 


D. EXTENSIONS TO MULTIPLE INTEGRALS 
The zero variance estimators above were derived in [2] 
only for single integrals. While the method can be extended 


to multiple integrals by nesting, nothing is gained since the 
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number of function evaluations increases exponentially with 
Bieenumber OL variables of integration. In this section, 
double integral extensions to the quadratic, cubic, quartic, 
and quintic zero-variance estimators, and triple integral 
extensions to the quadratic and cubic estimators are proposed, 
and the method of derivation is presented. Although not 
thoroughly tested on a wide range of integrands, it is be- 
lieved that this method of extension can, for multiple inte- 
grals with many variables of integration, provide an accurate 
Means of approximating the integral with fewer evaluations of 
the integrand than for the quadrature sum routines. 
aya 
For the double integral | | f(x,y) dxdy, consider the 
OO 
difference function which is a surface passing through the 
mmmereorners £(0, 0), £(0, 1), £(1, 0), and £(1, 1). This 


function is 
mean eee ety en lx) (lay) £070) + (1-x) yf (0,1) 
pee ay er ilo) + xyt( 1,1) (67) 
Integrating this difference function results in 


le 
[ [° tecesy-ee Cay) iaxay = £00481. 0/4600, 4 800) 
oO - 0 


(6 8) 


lieve: .*= 1S defined similar” to (55) as 


ie ree) 


I* ce 
Py (x) py (y) 


(62) 


and Py (X) and Py (y) HavemEpneommermeo: (59) solving for 
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fa, 7) 2467) and using (69) as the estimate for 


ie 7 1 Aa 
| [= £*(x,y) axay, the estimate to | | £(x,y)dxdy becomes 


on 0 ‘Ques 

i E009 +£(0,1L)+£ (1 fete 
Be] e(x,yrdxay = £10,0)48(0,1)+5(2,0)46(2,1) 

Oe - 0 


fp) x) ly) £00) (1) yt (0, 1) x(i-y) £01 ,0) ~xyf(1, 1) 
36xy (1l-x) (l-y) 


(70) 
which is the double integral extension to (56). Replacing x 
Paeeaei-x and y with I-y in (70) and averaging the resulting 
expression with (70) gives the double integral extension to 


mime cublLec Zero Variance estimator (62) 


Te ah 
| | nee = ees St) 0) +t, 1) 
0 -0 


f(x,y) +£(1-x,1l-y) -(1-x-y+t2xy) [£(0,0)+£(1,1) ] 
72xy (1L-x) (1-y) 


qs ieee Bee ae (71) 
72xy (1-x) (1-y) 


+ 


oa 


For higher order extensions, a difference surface is 
defined such that it passes through the a" end points of 
£(x'), where x Mca weve cton. Enis detference curve; which 
Toeaoolynomial in Gach ©f its variables, can be integrated 
Mieiecieatiy sto ugive a first approximation to the desired 
iitegcal. The estimate is then refined by sampling the f* 
function weighted by the distribution of the variables, as in 
(ope eunie the thiple extension to the quadratic estimator 


(56) becomes 
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ier pl - 
| | | £(x,y,2) axdyde = £(0,0,0)4£(0,0,1)+£(0, 1,0) +£(0, 1,1) 


eee ¢ 0 
, £(1,0,0)+£(1,0,1)+£(1,1,0)+f(1,1,1) | fn nS) 
8 216xyz (1-x) (1-y) (1-z) 


ae = | JES) (Ibs) CIL=2) £(0,0,0) -(1-x) i) Zi 0, Oe 
ZEGSyen =<) (oy) (So) 


ie mee yl 2) t 00,1, 0) 
216xyz(1-x) (l-y) (1-z) 


216xyz (1-x) (1-y) (1-z) 


aoa Zee jl — sy Zr { 1 1) 


* —“di6xyz(1-x) (I=y) (1-2) | (72) 


and the triple extension to the cubic estimator (62) becomes 


ier pl 
| | | ete advan ge 2 11) +6 (0,1 0) 48001) 
0 70 


0 


Meee 0,0) tf (1,0,1)+t(1,1,0)+£(1,1,/1) f(x,y,z)+f£(1-x,1-y,1-z) 


8 432xyz (1-x) (l-y) (1-z) 


(7 eee eee OOO (be 1) | 


ss 432xyz(1-x) (1-y) (1-2) 


i = (Zax e=VZEXy e070.) +E, 0) | 
432xyz (1-x) (1l-y) (1-z) 


. ceo Xa Zy tee COO) te le O74) | 
432xyz (1-x) (1l-y) (1-z) 


= (Ko K Za ty Za no) Oy LL) | 


: Ae oy ee le) 


(73) 


Bor the quartic and Gquincie double integral extensions, 
the difference surface must additionally pass through the 
four points which are peripherai midpoints and the internal 
Heep olnc Of the integrand. The double integral extension 


fomeeie quartic: Zero Vdariancegesctimator becomes 
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f(x,y) axdy..= 


[_ ie Or Oe 0, L) titi, Ore ( 1,1) 
0 36 


i iL at iL ibaa 
£(0,5) +£(1,5) +£(5,0) +£ (5,1) +4 (5,5) 
+ gps 4$__< e 
36 


Gee jeeee) a2 tly) (i=2y) £00, 0) 
2 


+ 
900xy (1-x) (1-y) (1-2x) * (1-2y) 


+4y (1-y) £(0,5) -y (1-2y) £(0,1)] 


900xy (1-x) (1-y) (1-2x) 7 (1-2y) ? 
D 2 
FOr ayaa ye) ele) LaDy) 


x(1-2x) [ (1-y) (1-2y) £(1,0) +4y (1-y) £(1,5) -y (1-2y) £(1,1)1 
— ey oan: aa nee an, 


r (74) 


900xy (1-x) (1-y) (1-2x) 7 (1-2y) 
Mme the double integral extension for the quintic zero 


Varlance estimator becomes 


fee L 
eee Oe Onder (nO el. 1) 
: I, £(x,y)dxdy = a, ee a 
£ (0,5) +£(1,5) +£ (5,0) 4£(5, L448 (5,5) 


+ 4[ : 
6 


meena (lx, 1—-y)i (1-2) (l-2y) (l-x-yr2xy) [£(0,0)+f(1,1)) 
2 





1800xy (1-x) (1-y) (1-2x) * (1-2y) 


(1-2x) (1-2y) (xty-2xy) [£(1,0)+£(0,1)] 


o D 


1800xy (1-x) (1-y) (1-2x) * (1-2y) 


~4y (1-2x) (1-y) [£(0,5) -£(1,5)] 
4- 





Hegel) (lay Glee <)> (1-2) 


~4x (1x) (1-2y) L£($,0) -£ (Fr 1) 1-32xy (1-x) (L-y) £(5/5) 


see a  -} 


ENO) Clon) (dl oos<) “lees 


Nie 
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For higher order extensions to the quartic and quintic 
estimators, the difference function must be defined so as to 


pass through the De corners, the ee peripheral midpoints, 


and the central midpoint of f(x? 


), where ce 1s a k vector. 
These extensions become difficult to derive er co Ene ange 
momoer Of points involved. 

The triple integral extension to the cubic estimator, (73) 
has been applied to the example of Chapter II. The results 
are contained in Chapter IV for comparison, and a printout 


of the program is contained in Appendix B. 
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EVs MeONCLUSIONS 


Chapter II presents the mathematics for obtaining a 
performance measure‘s probability density function for ex- 
plicit and implicit performance meaSures. Two integration 
schemes uSing Gauss quadrature sums were then discussed. 
Chapter III illustrated a crude Monte Carlo method for ob- 
taining a performance meaSure's distribution and then dis- 
cussed the confidence of such schemes. Monte Carlo schemes 
for evaluating single integrals were discussed. ‘Double and 
triple integral extensions to these methods were then derived. 

From obServation of the problem, it seems intuitive that 
feo oyorem S peretormance measures will be of the implicit 
type, and in general the performance meaSure will be a func- 
tion of several variables. In the process of evaluating the 
performance measure's distribution, the integral (1) must be 
evaluated. Evaluation by the quadrature sum methods, (25) 
and (29) require some a evaluations of the integrand, where 
n is the number of points in the quadrature sum, and k is the 
number of variables. This number can become quite large, and 
Since the computation time for such routines is roughly pro- 
portional to the number of times the integrand is evaluated, 
(lol es eee could take a good deal of time. However, this 
method iS quite accurate, and it is easily programmed for a 


large number of variables. 
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The reduced variance estimators reduce the number of 
evaluations of the integrand required. For a single esti- 
mate, the quintic extension requires evaluations at the os 
corners, the ees Pebtenerakemidooints an@*the 3 internal 
meres, Or a total of 3 +$(k+2)2*72 points, where again k 1s 
the number of variables. For additional estimates, only two 
additional points per estimate are required. Thus M estimates 
of the integrand would require 2M + 1 + VD De evaluations 
of the integrand. This can, 1n most cases be much less than 
that required for the quadrature sums. 

The drawback to the reduced variance scheme is that the 
geometric form of the integrand must be examined to insure 
@iatetne Major portion Of the surface will not be neglected in 
the sampling process. As illustrated in Figure 3, the dif- 
ference curve does not provide a good representation of the 
meron, wille splitting the curve into two Sections as in 
Figure 4 enables one to obtain a good representation of the 
function by the use of two difference curves. 

The generation of the random numbers with the probability 
denistty function AO Clas) pe) can be bothersome to program 
PomeCanmene Guintitc estimaror extension. The cubic estimator 
is much more easy to extend, and the density function 6x(1-x) 
1s more readily programmable. However more points are re- 
quired by the cubic to ensure as good accuracy as the quintic 
extension. 


Properly used, the Monte Carlo method of integral estima- 


tion can provide a more rapid method of determining a 
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Figure 4. 
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performance measure's density function with reasonable ac- 
curacy than other methods which require more standard inte- 


gration routines. 


Pee COMMENDATIONS FOR FURTHER STUDY 

Monte Carlo schemes seem to be attractive for reducing 
the amount of time required in the evaluation of the multi- 
ple integrals required in evaluation of a performance 
measure's distribution. Higher order zero variance estima- 
tors could be generated by fitting a difference curve to the 
points which are the zero's of a Legendre volynomial and 
evaluating the required probability density function to 
achieve a zero variance estimator for integrands of order 
2n-l, where n is the order of the Legendre polynomial. Other 
schemes of generating geist MUNDessSwOtepatelcular aust ribu- 


tions could as well be investigated. 


PPeeekaoUles OF APPLICATIONS TO EXAMPLE 

Hieweumibati ve O1Strlbputelon function of the upper 3 db 
frequency of the series R-L-C circuit of Figure 1 has been 
obtained by four methods. 

First, using the analytical expressions (11) and (12) for 
the integral (6), the probability density function, equation 
(15) was obtained by evaluating the integral with a 6-point 
Gauss Hermite quadrature sum routine. The cumulative dis- 
tribution function was then obtained from the probability 
density function by rectangular integration. The resulting 


distribution has been plotted in Figure 5, and will serve as 
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the reference since integration of explicit functions is 
generally more accurate than integration of implicit 
muMectLons. 

Using the same performance measure as for the explicit 
fase, but Seen eing that the value of the performance meas- 
ure was to be obtained implicitly, integral (18) was evalu- 
ated uSing 6-point Gauss Legendre quadrature sums. As can be 
seen in the plot of the cumulative distribution function for 
this case, Figure 6, the distribution is very nearly the same 
as that obtained in the explicit case. 

Then for the crude Monte Carlo case, 1000 samples of the 
performance measure (10) were taken with the values of the 
three variables being obtained from a Gaussian random number 
generator. As can be seen in Figure 7, the resulting cumula- 
tive distribution function 1S very nearly the same as that 
obtained in the explicit case. 

The extended cubic zero variance estimator equation (73) 
was then applied to the integral (18). The required distri- 
butions of the components were generated as explained in 
Appendix A. Twenty approximations for each value of the per- 
formance measure averaged to obtain the cumulative distribu- 
tion plotted in Figure 8. Two discrepancies are observed 
here. First, the cumulative distribution exceeds the maximum 
value of 1, and secondly, there is an apparent discontinuity 
in the curve. Since the entire curve has a higher value than 
that of Figure 5, it 1S assumed that the difference surface 


estimation was too large. 
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Although the example did not prove the worth of the 
Monte Carlo method, a count of the required number of evalu- 
ations for high multiple integrands for a 6-point Gauss 
Legendre routine as compared to that required in the extended 
cubic estimator averaging 100 approximations for each value 
of the performance measure will illustrate its merits. For 
a multiple integral with 10 variables of integration, the 
6-point Gauss Legendre method would require eee Beis Tae 
of the integrand. For the same integrand, the cubic exten- 

10 


Sion would require 200 + 2 evaluations. The savings in 


computation time should be apparent. 
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Appendix A 
RANDOM NUMBER GENERATION 


A. UNIFORM RANDOM NUMBERS 

A truly uniform random number cannot be generated on a 
digital computer. However, methods such as the power residue 
[4, 10] method can generate a sequence of mumbers which have 
the properties of uniformly distributed random numbers. 

The method of power residues 1s based on modular arith- 
metic. If Xs is an arbitrary integer in the sequence of 
pseudo-random numbers to be generated, the next number gener- 


ated is determined by 


ei laa X.p (modulo M) =A) 
where p 1S any prime integer such that oy Momeere 16e tat iicmete, 
Paeemodular base of the computer's arithmetic unit. The 
number xX. /M 1s then an approximation to a random variable 
Peeweadeune form GiStribution om the interval (0,1). Utilizing 
peoereatea! scolmpucer wrth va 32 bit register (31 number bits 
Psomodems tonite), a valle Of p Of 65,539 for an M of 
2,147,483,64/ = Bora a sequence of ae numbers can be gener- 
ated before any number of the sequence is repeated [4]. 

Since the digital computer automatically performs modular 
arithmetic on integers, with the modulus being determined by 


the size of the registers in the arithmetic unit, uniform 


pseudo random numbers can be generated quite rapidly. 
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B. GAUSSIAN RANDOM NUMBERS 

GauSSian random numbers are most commonly generated uti- 
lizing uniformly distributed pseudo-random numbers. If 
Ry, 1=1,2,...,n are independent samples from a density func- 
tion uniformly distributed over the interval Uo a then a 


normally distributed random variable with 0 mean and vari- 


ance of 12/n can be approximated by [3, 4, 10] 


— 12 , 
=) CR) = 6 (A-2) 
LIL 
meena — 12, Z4Z will have 0 mean and variance of l. 


eee ox({!—-x) DISTRIBUTION 

lea R. is an independent sample from a random variable 
Pmbeonrmiy Qustriputea on the interval (0,1), for 1 = 1,2,3, 
and if the three samples are rearranged in order of magni- 


tude such that 


Re RoR 


1 5 (A-3) 


3 


then Ry fismiomeodomidiy ensituy function 6x({l-x) [3]. 


D. 30x(1-x) (1-2x)* DISTRIBUTION 


te Rie Ros Ro, Rae and R. are independent samples froma 


s) 
random variable uniformly distributed over the interval (0,1), 


and are arranged so that 


|< [R, : =| See | Reo =| (A-4) 


Nib 


[Ry = 


and if S is chosen such that S = Ry with probability 3/4 and 


S = R, with probability 1/4, then S has probability density 
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puncte Lon 30x (1-x) (1-2x)? tweets Gensity function, it was 


discovered, requires some complex programming to generate. 


eee) OTHER DESTRIBULTIONS 
Let Y be a uniformly distributed random variable over the 


mmtcerval (0,1), and let 


Ss 
y = F(x) = | lov ie) she (A-5) 
wmente Mit) = L for 0 t 1, and h(t) = 0 otherwise. Since 
Xx 
F(x) = | Py (x) (A-6) 


—OO 


where Py (x) MomeicewoTOodeLIaty Gensity function for X, if 
Hix) Or Py (x), the desired function, are known explicitly, 


jmen the inverse function 
ues oF {y) Cae) 


will generate x according to its desired density function 

Py (x) Pieweeiocst often Py (x) wWoll be- a polynomial fitted to 
some statistical data which represents the distribution of 
the variable X. Since this polynomial is normally of a high 
order, some difficulty can result from this method in that 


EiemgZeis@cmo: (A=G) May be difficult to find. 
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APPENDIX 86 
COMPUTER PROGRAMS 
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